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A precessing jet model for the PN K 3-35: simulated 
radio-continuum emission 
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ABSTRACT 

The bipolar morphology of the planetary nebula (PN) K 3-35 observed in radio- 
continuum images was modelled with 3D hydrodynamic simulations with the adap- 
tive grid code YGUAZU-A. We find that the observed morphology of this PN can 
be reproduced considering a precessing jet evolving in a dense AGB circumstellar 
medium, given by a mass loss rate M csrn — 5 x 10 -5 M0yr -1 and a terminal veloc- 
ity v w = 10 kms -1 . Synthetic thermal radio-continuum maps were generated from 
numerical results for several frequencies. Comparing the maps and the total fluxes 
obtained from the simulations with the observational results, we find that a model of 
precessing dense jets, where each jet injects material into the surrounding CSM at a 
rate Mj = 2.8 x 10 -4 M yr _1 (equivalent to a density of 8 x 10 4 cm -3 ), a velocity 
of 1 500 km s -1 , a precession period of 100 yr, and a semi- aperture precession angle 
of 20° agrees well with the observations. 

Key words: hydrodynamics - methods: numerical - ISM: planetary nebulae: indi- 
vidual K 3-35 - ISM: jets and outflows - radiation mechanisms: thermal 



1 INTRODUCTION 

The evolution between the end of the asymptotic giant 
branch (AGB) and the planetary nebula (PN) phases has 
for a long time been a poorly understood link in the late 
stages of intermediate-mass stars (1-8 M©). It is in this 
transition phase where the fast stellar wind of the emerging 
PN interacts with the slow wind from its precursor AGB 
star (Kwok, Purton, & Fitzgerald 1978), shaping the final 
morphology of the PN. 

Other ingredients that can contribute to the final shape 
of a PN are the presence of interacting binary stars (e.g. 
Morris 1987; Balick & Frank 2002; Soker & Bisker 2006) or 
magnetic fields (e.g. Garcfa-Segura & Lopez 2000; Blackman 
et al. 2001; Garcfa-Segura, Lopez & Franco 2005). Thus, 
the study of objects in this transition phase can give us 
important clues about the physical mechanisms responsible 
for the different morphologies observed in PNe. 

K 3-35 is a very young PN with a characteristic S- 
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shaped emission morphology that suggests the presence of 
precessing bipolar jets (Aaquist & Kwok 1989; Aaquist 1993; 
Miranda et al. 2000; 2001). The detection of OH and H 2 
maser emission as well as the presence of CO and HCO + 
emission suggest that K 3-35 departed from the proto-PN 
phase only a few decades ago (Miranda et al. 2001; Tafoya 
et al. 2007). 

Miranda et al. (2001) estimate a dynamical age of ^ 
15 years for the ionised inner core, which expands at ^25 
km s _1 . For the jets, assuming a modest jet velocity of ^100 
km s _1 , an age of ^800 years is obtained. Therefore, the jet 
formation in K 3-35 occurred during the proto-PN phase. In 
this paper we show that the radio morphology of K 3-35 can 
be explained by a precessing jet model. 



2 THE PRECESSING JET MODEL 

The PN K 3-35 belongs to the point-symmetric PN 
group because it shows a bipolar structure, with an 'S' 
shaped morpho l ogy i n radio - continuum maps at 8.3 GHz 
(iMiranda et al] l200lL Il998l : lAaauistl Il993l : iGomez et 



2003). 
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This kind of morphology is more often observed in 
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Figure 1. Scheme of the initial conditions for the simulation. 
The xyz— frame corresponds to the computational domain. The 
physical system is given by the x'y'z'— frame, which is obtained 
by rotating the yz-plane around x-axis by 36° with respect to the 
plane of the sky (xz-plane). 



galactic and extragalactic jets. Several theoretical studies 
have employed precessing jets for modelling these objects 
and explaining the 'S' shape observed at different frequen- 
cies. 

F or instance, in the galactic case, iMasciadri et al.l 

(2002) carried out 3D numerical simulations of the HH 34 
jets, and compared jet models wi th and without precession. 
In the field of extragalactic jets, iRodriguez-Martmez et al.l 
(2006) made a 3D hydrodynamical simulation of a precess- 
ing jet moving into an intragalactic cluster medium, and 
compared their results with the Abell 1795 galaxy cluster. 
The similarity between these jets and point-symmetric PNe 
has recently motivated models of this type of PNe as jets, 
considering precession, density and/or variable ejection ve- 
locity, etc., in order to explain their kinematics, emission, 
and morphology. An exce llent review of these works is given 
bv lSoker Biskerl (|200d ). 

As an example, we can mention the observational 
study of th e prot o-planetary nebula (PPN) Hen 3-1475 by 
iRiera et al.l (|2QQ3h . In Ha, [Nil], and [OIII] imag es, this PPN 
exhibi ts several knots forming an 'S' shape. Also lRiera et al .1 

(2003) made radial velocity and proper motion studies, 
showing that the knotty structure of Hen 3-1475 can be ex- 
plained by means of a jet with a variable ejection velocity 
and p recession. Based on these hypotheses, Velazquez et al.l 

(2004) have computed 2D (axisymmetric) and 3D hydrody- 
namical simulations of such a jet, showing that this model 
can reproduce both the observed morphology and the kine- 
matical characteristics of the outflows of Hen 3-1475. 

In the present paper, we model the PN K 3-35 with 3D 
numerical simulations of a precessing jet. Due to the fact 
that the bipolar structure of K 3-35 looks almost continuous, 
we do not find necessary to consider also a time variability 
for the jet velocity and/or density. 

In PNe, the surrounding circumstellar medium (CSM) 
is strongly modified by the central star, which injects mate- 
ri al and inc r eases the density of the CSM. Based on the work 
of lMellemal(|l995h . this perturbed CSM can be described by: 



p(R) po 



1 -exp(-2/3sin 2 6) 
1 - exp(-2/3) 



(Ro/R) 2 



(1) 



where R — x 2 + y 2 + z 2 is the distance from the jet source 
on a 3D Cartesian frame. The angle is measured with re- 
spect to the equator, which we have chosen to coincide with 
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Figure 2. Temporal evolution of the column density at integra- 
tion times t = 20, 30, 40, 50 yrs. These maps were obtained inte- 
grating the density along the x-axis. At t=40 yr, the jet achieves 
an extension of 1.8 x 10 17 cm, which is equivalent to 2.4", i. e., 
the angular extension of bipolar structure of PN K 3-35 (for an 
assumed distance of 5 kpc). Both axes are labelled in units of 
10 17 cm, and the bar on the right represents the logarithmic gray- 
scale, given in units of cm -2 . 



the x2/-plane, yielding = tan -1 ( \/ x 2 + y 2 / ' z). Eq. [T] de- 
scribes the stratification of the surrounding CSM produced 
by a star in AGB phase. The parameter 5 in Eq. [1] is re- 
lated to the ratio of the density at the equator (xy-plane) 
and at the pole (z-axis), and (3 determin es how the den - 
sity varies from the equator to the pole (|Mellemal [T995). 
The value of po is calculated from the mass loss rate as : 
po = M/(47r Rq Venv), where v env is the expansion velocity 
of the AGB remnant. 



3 THE NUMERICAL SIMULATIONS 

3.1 Initial numerical simulation setup 

The numerical simulations were car ried out employin g the 
3D adaptive grid code yguazu-A (Ra ga et all l2000h . To- 
gether with the gas-dynamic equations, several rate equa- 
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10 cm along the 



direction, respectively. These high 



densities are supported by the recent detection of HCO + 
toward K 3-3 5, that reveals the presence of high density 
molecular gas (jTafova et alJl2Q07h . 



tions for atomic/ionic species (H i-ii, He i-iii, C II- IV, and O 
i-iv) are integrated in order to calculate a non-equilibrium 
ionisation cooling function for the gas (for temperatures 
above ~ 10 5 K the cooling rate is approximated by an 
equilibrium ionisation parameterised function). The gas- 
dynamic equations are integr a ted w ith the "flux- vector split- 
ting" algorithm of I Van Leerl (|l982l ). Details of the reaction 
and cooling rates included in the c ode are described in the £ 
appendix of iRaga et~ai1 (|2QQ2h and iRaga et al.1 (|20Q7h . ^ 

The simulations were performed in a five level, binary 2 
adaptive grid, with a maximum resolution of 5.9 x 10 14 cm x _ 

(along the three axes). The computational domain had an 
extent of 1.5 x 10 17 cm in the x- and y- direct ions, and of 
3 x 10 17 cm in the ^-direction. 

Both the cone precession axis and the CSM density dis- 
tribution were rotated around the x-axis, by 36°with respect 
to the plane of the sky (?/£-plane, see Fig [T]) . This inclination 
angle was found bv lUscanga et all ([2007h , who modelled the 
H2O maser emission of the accretion disk around the K 3-35 
central star. The configuration is shown in Fig. [1] 

The jet is placed at the centre of the computational do- 
main with an initial radius r- } — 1.5 x 10 16 cm (equivalent 
to 30 pixels at the maximum grid resolution) and a length 
Zj = 6 x 10 15 cm. Initially (integration time t — yr), the 
jet axis lying on the x'z'-plane and starts to move (coun- 
terclockwise) into a precession cone, with a semi-opening 
angle a — 20°. At t — yr, the jet velocity components are 
Vj = (464., 0., 1426.) km s" 1 , in the x'y'z'— frame (see fig. 
□J, or vj = (464., -764., 1204.) km s" 1 , in the xyz- frame 
(the computational domain frame). After several tests, a pre- 
cession period of r p = 100 yr was chosen. Two runs were 
computed considering that the jet material is injected as 
singly ionised or neutral gas. K 3-35 has a bright centre 
when it is observed at radiofrequencies, because the cen- 
tral source strongly ionises the CSM gas. This process is not 
taken into account in our simulations. Since we continuously 
impose the jet inflow conditions within a cylindrical volume 
(see §2), when we calculate synthetic radio-continuum maps 
we obtain a strong central emission component (coming from 
this cylinder) for the initially ionised jet simulation, giving a 
good morphological agreement with the observations. How- 
ever, this central component in the simulated radio emission 
is exclusively a result of the chosen jet inflow conditions, 
and for this reason only synthetic maps obtained for the 
run with the "neutral jet" will be employed for calculating 
total fluxes. 

A good fit with the observations is obtained when 
the velocity and mass loss rate of the jet were set to 
1,500 km s" 1 and Mj = 2.8 x 10" 4 M© yr" 1 , respec- 
tively. These parameters let an equivalent density pj of 
8 x 10 4 cm" 3 . The neutral CSM is modelled by Eq[Q with 
a mass lose rate M — 5 x 10" 5 Mq yr" 1 , v w — 10 6 cm s" 1 , 
and the parameters S and /3, were set to 0.1 and 3, respec- 
tively. These parameters give an almost isotropic density 
distribution with a contrast of 20% between the density at 
the equator and the poles. From this relationship, densities 
between 7.5 x 10" 19 g cm" 3 and 1.7 x 10" 20 g cm" 3 are ob- 
tained for distances from the source of 1.5 x 10 16 cm and 
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Figure 3. Simulated thermal radio-continuum maps obtained for 
a frequency of 8.3 GHz. The left panel is the emission map gener- 
ated by integrating the emission coefficient along the y-axis, the 
right panel is the result of integrating along the x-axis. The bar 
on the right shows the logarithm gray-scale of the radio thermal 
emission in units of erg s _1 cm -2 sr -1 . Both spatial axis are in 
units of 10 17 cm. Strong emission is observed in the centre due 
to the fact that this is the region where a singly ionised jet inflow 
condition is continuously imposed. 



Taking into account the CSM density distribution, and 
the jet density and velocity, it is possible to estimate the 
speed of the leading 'working surface' (i. e. the head) of the 
jet : 



VWpcsm (i?,6>) 

VWpcsm(fl,0) + l' 



(2) 



where v$ is the jet injection velocity, and pj and p csm are 
the jet and CSM densities, respectively. It is straightforward 
to show from Eqs. J2} and JT]) that for large distances R 
from the jet source, v ws increases. For instance, considering 
the values for pj, and p CS m, speeds of 4.4 x 10 2 km s" 1 
and 1.1 x 10 3 km s" 1 , are found at R = 5 x 10 16 cm and 
R = 10 17 cm , respectively. 



3.2 Calculation of the thermal radio-continuum 
emission 

To generate synthetic thermal radio-continuum maps, the 
radiative transfer equation was integrated along different 
lines of sight. 

First, the optical depth r u (l) was calculated by inte- 
grating the absorption coefficient k, u along lines of sight : 
r u (l) — j Q K v dl (with / being the spatial coordinate along 
the line of sight). The absorption coefficient k u is given by 
(jRvbicki fe LiditmaiJll97& IOsterbrocklll989h : 



0.018 T 



-3/2 



n e mv g s , (3) 

where # ff = 0.5513 (/n(T 3/2 z/" 1 ) + 17.7) is the Gaunt factor. 
We obtain directly from our simulations T, n e , and m (the 
spatial distributions of temperature, electronic, and ionic 
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21.7 


13.4 
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24.0±2.4 
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19.1 
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Table 1. Total flux densities (in units of mJy) at different fre- 
quencies, obtained from the simulated images at t = 40 yr (second 
column), and at t = 50 yr (third column). In this calculation, a 
distance of 5 kpc has been assumed to PN K 3-35. These flux 
densities have been obtained from the simulation where the jet 
material is imposed as neutral gas, in order to avoid an unrealis- 
tic strong emission from the centre. For a direct comparison with 
observation s, the fourth c olumn gives the "halo" observed fluxes, 
taken from lAaauistl (|l993h . 

number densities, respectively), which are used in Eq. (|3} to 
calculate k v (in units of cm -1 ). 

Then, the intensity as a function of frequency can be 
obtained as : 

Iu= f ju(n ei n tl T) e~ TAl) dl, (4) 
Jo 

where j u is the emission coefficient (which is a func- 
tion of T, n e , and m). The intensity l v has units of 
erg s -1 sr -1 cm -2 Hz -1 . Numerically, the emission coef- 
ficient can be calculated using the follow ing equation 
(|Rvbicki fe Light man|[l979l : IOsterbrocklll989h : 

c /iiio w in-39 rp-1/2 -hvlkT / K \ 

> — 5.4113 x 10 n e m T 1 e 1 g^, (5) 

where h and k are the Planck's and Boltzmann's con- 
stants, respectively. The coefficient j u has units of 
erg s -1 sr -1 cm -3 Hz -1 . 

Synthetic maps at different frequencies were generated 
combining the results from the 3D simulation with Eqs. ((3}, 
Q, and (0. The angularly integrated fluxes (in mJy) corre- 
sponding to different frequencies were then determined from 
the si mulated maps, assuming a distance to the source of 
5 kpc (iZhandll995h . 



4 RESULTS 

With the physical conditions described in §3.1, a 3D numer- 
ical simulation was carried out. Figure [2] shows the tempo- 
ral evolution of the column density obtained for integration 
times lying in the range [20,50] yrs. The stratification of the 
surrounding CSM is clearly observed in all maps. 

Filamentary structures are observed at the head jets, 
due to precession movement. At different times, jet material 
is injected at different angles. These jet material sweep up 
the CSM which have a density stratification as function of 
the angle (see Eq. [TJ), generating the bright filaments ob- 
served in the column density maps. 

At a t — 40 yr integration time, the jet reaches a total 
length of 1.8 x 10 17 cm. This size is equivalent to 2.4" (i. e., 
the angular exte nsion of K 3-35 on th e plane of the sky) , at a 
distance of 5 kpc. lMiranda et all (|200lh have made estimates 



of the dynamical ages of both the inner ionised core and the 
extended emission (jets) of the PN K 3-35. Based on an 
expansion rate of ~ 25 km s _1 taken from the He II lines, 
they obtained a dynamical age of ^ 15 yr for the core, while 
assuming a jet velocity of ~ 100 km s _1 , a dynamical age of 
~ 800 yr is deduced for the jets. Our simulations required a 
much larger jet velocity to explain the radio emission of the 
lobes in the PN K 3-35, thus yielding a short dynamical age 
for the jets, which is comparable with that of the core. This 
is interesting because it may suggests that the process that 
originated the ionisation of the jets and the core occurred 
almost simultaneously. 

Synthetic maps were generated (as described in section 
3.2) for the following frequencies: 1.5, 2.4, 5, 8.3, 10, 15, 
and 22 GHz. Very similar morphologies are obtained for all 
of these frequencies. Figure [3] shows the maps obtained at 
8.3 GHz. The left panel of Figure [3] displays the simulated 
radio- continuum emission when the j's.sghz (Eq. [5} is inte- 
grated along the ?/-axis, giving the emission on the xz-plane. 
The right panel of Figure [3] shows the 8.3 GHz emission on 
the 2/z-plane (i. e., the integration of j8.3GHz was carried 
out along the x-axis). The map on the right was shown to 
illustrate how the morphology strongly depends on the pro- 
jection chosen. In the xz- project ion, the radio emission has 
an 'S' shape (which is the morphology of K 3-35), while 
in the yz-projection, we have a 'Z' shape. These maps are 
displayed in units of erg cm -2 s _1 sr -1 . 

Comparing the left panel of Figure [3] with the column 
density map at t = 40 yr of Figure [2] we see that the radio- 
continuum emission traces very well the precessing jet mor- 
phology. 

It is also possible to compare simulated and observed 
total fluxes. However, since we are not simulating the radio 
emission from the central region, we can only compare the 
radio e mission fr o m the lobes or jets, which were called the 
'halo' in lAaquistl ([l993). From each simulated radio map we 
measured the corresponding total flux from these structures. 
These fluxes are listed in Table 2] where one can see that the 
total fluxes in the simulations (at all frequencies) are of the 
order of the observed ones. At an integration time of 40 yr, 
the simulated flux gives an optically thin behaviour for the 
spectrum. On the other hand, for integration times of 50 yr, 
the simulated flux has a maximum at a frequency between 
5 and 8.3 GHz, in good agreement with the Aaquist (1993) 
results. At frequencies higher than 8.3 GHz, the spectral in- 
dex is of —0.09, which corresponds to optically-thin thermal 
free- free emission. 

A direct comparison between both observed and simu- 
lated 8.3 GHz maps is shown in Figure [4] The simulated map 
was convolved with a synthesised beam of 0.2" x0.2 7/ (similar 
to the beam in the observations), and scaled to yield units 
of mJy beam -1 . The same scale for both images have been 
employed. Also, the synthetic map was tilted in 48° (with 
respect to z— axis, see Fig. 1) in order to do the position 
angle (P.A) of the jet direction (in the xz— plane) equal to 
65°, the observed P.A. ([Miranda et al.|[2QQlh . 



5 SUMMARY 

Several 3D hydrodynamical simulations were carried out em- 
ploying the adaptive grid code yguazu-A, in order to model 
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both the morphology and the thermal radio emission of the 
planetary nebula K 3-35. 

After analysing our results, we found that the bipolar 
structure of this PN can be described as the result of the 
interaction of a dense jet (with an initial number density of 
8 ~ 10 4 cm -3 or M 3 = 2.8 x lO~ 4 M 2/r _1 ) moving into a 
dense environment, previously swept up by the AGB wind 
of the central star. 

The 'S' morph ology shown by K 3-35 in 8 . 3 GHz radio- 
conti nuum images ([Miranda et al.| [l998. 2Q01; [G6mez et ah! 
2003) can be reproduced if the modelled jet precesses with 
a period of 100 yrs on a cone with a half-opening angle of 
20°. 

For an integration time of 40 yrs, the simulated jet 
has a total length of 1.8 x 10 17 cm, which i s equivalent to 
2.4 // (lMiranda et al.|[200ll : iGomez et alJl2003h considering an 
estimated distance of 5 kpc to K 3-35 and also the orienta- 
tion of this object ([Uscanga et aUkOOTlV This time is almost 
20 times smaller than the one given bv lMiranda et al.l ([200 lh 
for the jet, where a slower velocity was assumed. However, 
it is only 2.7 time s larger than the dy namical age of the in- 
ner core (also in iMiranda et all l200lh , suggesting that the 
two events are more related than previously thought, and 
further supporting the idea that K 3-35 is a young object. 

Synthetic radio- continuum maps were generated from 
our numerical results. These maps show that the predicted 
morphologies and fluxes are in reasonable agreement with 
the observations. At an integration time of 40 yrs, the ob- 
tained spectral index is the one of optically thin emission. 
For an integration time of 50 yrs, the observed change of the 
spectral index with frequency (Aaquist 1993) is also repro- 
duced by our simulation. A direct comparison between the 
observational and numerical results is given in Fig. [4] where 
we show the observed and synthetic (for an integration time 
of 40 yrs) 8.3 GHz radio maps. 

We must note that the values for the velocity and mass 
loss rate employed in the simulation for the jet and the CSM 
seem to be rather high. It is difficult for AGB and post- AGB 
starts to launch jets at velocities of the order of 1 500 km s _1 
(although Riera et al. 2003 have reported this kind of veloc- 
ity for the outflows of the PPN 3-1475). Besides, the mass 
injection is also a bit extreme, in 40 years both jets have 
injected 0.02 M into the surrounding CSM. This scenario 
might be explained in terms of a binary system, if the jet 
is produced by a companion accreting material (at a rate 
about ten times higher than Mj) from a massive star (the 
AGB progenitor that produced the density distribution of 
the circumstellar material) . The primary in the last 40 years 
has lost 0.2 M©. This means that it started with an enve- 
lope containing this mass, which appears to correspond to an 
AGB star rather than a post- AGB star (even though at the 
present time the star has already evolved to the planetary 
nebula stage). 

Furthermore, the CSM is very dense. With the param- 
eters employed, in a radius of 10 17 cm, a mass of 1 M© is 
contained, implying a massive PN. There is observational ev- 
idence that favours a very dense surrounding CSM. However, 
the total mass derived from HCO + observations (Tafoya et 
al. 2007) is quite low (~ 0.017 M ), so that the molec- 
ular emission from this molecule would be confined to a 
small, possibly shock excited volume, in order to be consis- 
tent with our much higher mass CSM. These values imply 



that this scenario would be plausible if it is a short lived 
event. Clearly, a better determination of these parameters 
could be done with proper motion studies and better dis- 
tance determinations. 

Notwithstanding the extreme values for the employed 
parameters, it is important to note that a simple model 
of a precessing dense jet moving into an also dense CSM, 
successfully reproduces the observed morphology of the PN 
K 3-35, obtaining a qualitatively and quantitatively good 
agreement between the model predictions and the observa- 
tions, although this scenario would be a short lived event. 
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Figure 4. Comparison between observational and numerical re- 
sults. The top panel shows the 8.3 GHz radio co ntinuum image 
of K 3-35 ([Miranda et al.ll200ll : [Gomez et al.ll2003l ). The bottom 
panel shows the tilted synthetic radio emission map at the same 
frequency, obtained for an integration time of 40 yr. An inclina- 
tion of 36° between the precession axis and the plane of the sky 
was considered. The simulated map was smoothed with the ob- 
served beam of 0.2 /; x 0.2 /; and it is shown in the same scale that 
the observed one. The contours corresponds to the levels 0.18, 
0.31, 0.56, 1.0, 1.78, 3.16, 5.62, and 10.00 mJy beam" 1 . The ar- 
rows indicate the jet direction (projected on the plane of the sky) 
for integrations times of and 40 yr. The last o ne is coincident 
with the P. A. measured bv lMiranda et al.l (|200lh . 
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